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We study the magnetic braking and viscous damping of differential rotation in incompressible, 
uniform density stars in general relativity. Differentially rotating stars can support significantly 
more mass in equilibrium than nonrotating or uniformly rotating stars, according to general rela- 
tivity. The remnant of a binary neutron star merger or supernova core collapse may produce such a 
"hypermassive" neutron star. Although a hypermassive neutron star may be stable on a dynamical 
timescale, magnetic braking and viscous damping of differential rotation will ultimately alter the 
equilibrium structure, possibly leading to delayed catastrophic collapse. Here we treat the slow- 
rotation, weak-magnetic field limit in which E rot -C -Emag <§C W, where E Iot is the rotational kinetic 
energy, -E ma g is the magnetic energy, and W is the gravitational binding energy of the star. We 
assume the system to be axisymmetric and solve the MHD equations in both Newtonian gravitation 
and general relativity. For an initially uniform magnetic field parallel to the rotation axis in which 
we neglect viscosity, the Newtonian case can be solved analytically, but the other cases we consider 
require a numerical integration. Toroidal magnetic fields are generated whenever the angular veloc- 
ity varies along the initial poloidal field lines. We find that the toroidal fields and angular velocities 
oscillate independently along each poloidal field line, which enables us to transform the original 
2+1 equations into 1+1 form and solve them along each field line independently. The incoherent 
oscillations on different field lines stir up turbulent-like motion in tens of Alfven timescales ("phase 
mixing"). In the presence of viscosity, the stars eventually are driven to uniform rotation, with the 
energy contained in the initial differential rotation going into heat. Our evolution calculations serve 
as qualitative guides and benchmarks for future, more realistic MHD simulations in full 3+1 general 
relativity. 

PACS numbers: 04.40.Dg, 95.30.Qd, 97.60.Jd 



I. INTRODUCTION 

New-born neutron stars formed from core collapse su- 
pernovae, accretion induced collapse of white dwarfs, and 
coalescence of neutron star binaries are likely to be differ- 
entially rotating [1-6]. Differential rotation will cause a 
frozen-in poloidal magnetic field to wind up and generate 
a strong toroidal field. The back-reaction of the magnetic 
stresses on the fluid motion - magnetic braking - will de- 
stroy the initial differential rotation and can result in a 
significant change in the structure and dynamics of the 
star. 

Differentially rotating neutron stars can support sig- 
nificantly more rest mass than their nonrotating or uni- 
formly rotating counterparts, creating "hypermassive" 
neutron stars [16, 17]. Such hypermassive neutron stars 
can form from the coalescence of neutron star bina- 
ries [5, 6]. The stabilization arising from differential 
rotation, although expected to last for many dynami- 
cal timescales, will ultimately be destroyed by magnetic 
braking and/or viscosity [16, 38]. These processes drive 
the star to uniform rotation, which cannot support the 
excess mass, and can lead to "delayed" catastrophic col- 
lapse, possibly accompanied by some mass loss. 

Conservation of angular momentum implies that the 
nascent neutron stars formed from core collapse super- 
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novae or accretion induced collapse of white dwarfs are 
rapidly and differentially rotating. These neutron stars 
may develop a nonaxisymmetric bar-mode instability 
that might be detectable as a quasi-periodic gravitational 
wave by gravitational wave detectors, such as LIGO, 
VIRGO, GEO and TAMA. The dynamical bar-instability 
usually develops when the ratio /? = T/W is high enough, 
e.g. (3 > 0.25 - 0.27 [4, 18], where T is the rotational ki- 
netic energy and W is the gravitational binding energy 
(however, for extreme degrees of differential rotation a 
bar-instability may also develop for stars with very low 
f3 [19-21]; moreover, an m — 1 one-armed spiral mode 
may also become dynamically unstable for stars with a 
very soft equation of state and a high degree of differen- 
tial rotation [31, 32]). A hot, proto-neutron star formed 
a few milliseconds after core bounce may not have a suf- 
ficiently high value of j3 to trigger the dynamical bar- 
instability immediately. However, the instability might 
develop after about 20 seconds, after which neutrinos will 
have carried away most of the thermal energy, causing 
the star to contract further and (3 to exceed the thresh- 
old value for the bar-instability. If the proto-neutron 
star has a strong magnetic field (B > 10 12 G), magnetic 
braking could significantly change the angular momen- 
tum distribution during the neutrino cooling phase and 
might suppress the bar-instability [4]. On the other hand, 
the secular bar-instability could develop at much lower 
f3 [7] on a longer timescale. Once again, a small, frozen- 
in seed magnetic field could be wound up to sufficiently 
high strength by differential rotation to suppress the sec- 
ular instability over this timescale. 
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Short-duration gamma-ray bursts (GRBs) are thought 
to result from binary neutron star mergers [22], or tidal 
disruptions of neutron stars by black holes [2]. Long- 
duration GRBs likely result from the collapse of rotating, 
massive stars which form black holes, followed by super- 
nova explosions [15, 23]. In current scenarios, the burst 
is powered by the extraction of rotational energy from 
the neutron star or black hole, or from the remnant disk 
material formed around the black hole [24] . Strong mag- 
netic fields provide the likely mechanism for extracting 
this energy on the required timescale and driving colli- 
matcd GRB outflows in the form of relativistic jets [25- 
27] . Even if the initial magnetic fields are weak, they can 
be amplified to the required values by differential motions 
or dynamo action. 

The r-mode instability has recently been proposed as 
a possible mechanism for limiting the angular velocities 
in neutron stars and producing observable quasi-periodic 
gravitational waves [39-42]. However, preliminary calcu- 
lations (see Refs. [43-45] and references therein) suggest 
that if the stellar magnetic field is strong enough, r-mode 
oscillations will not occur. Even if the initial field is weak, 
fluid motions produced by these oscillations may amplify 
the magnetic field and eventually distort or suppress the 
r-modes altogether (r-mode oscillations may also be sup- 
pressed by nonlinear mode coupling [46-48] or hyperon 
bulk viscosity [49-51]). 

In a different context, supermassive stars (SMSs) may 
form in the early universe, and their catastrophic col- 
lapse may provide the seeds of supermassive black holes 
(SMBHs) observed in galaxies and quasars (see Refs. [28- 
30] for discussion and references). If an SMS maintains 
uniform rotation as it cools and contracts, it will ulti- 
mately arrive at the onset of a relativistic radial instabil- 
ity, triggering coherent dynamical collapse to an SMBH 
and giving rise to a burst of gravitational waves [33, 34]. 
If an SMS is differentially rotating, cooling and contrac- 
tion will instead lead to the unstable formation of bars 
or spiral arms prior to collapse and the production of 
quasi-periodic gravitational waves [35, 36]. Since mag- 
netic fields and turbulent viscosity provide the principal 
mechanisms that can damp differential rotation in such 
stars [37, 38], their role is therefore crucial in determining 
the fate of these objects. 

Motivated by the importance of magnetic fields in dif- 
ferentially rotating relativistic stars, Shapiro performed 
a simple, Newtonian, MHD calculation showing some of 
the effects of magnetic braking ([38], hereafter Paper I). 
In his simplified model, the star is idealized as a differ- 
entially rotating, infinite cylinder consisting of a homo- 
geneous, incompressible conducting gas. The magnetic 
field is taken to be purely radial initially and is allowed 
to evolve according to the ideal MHD equations (flux- 
freezing) . The calculation shows that differential rotation 
generates a toroidal magnetic field, which reacts back 
on the fluid flow. In the absence of viscous dissipation, 
the toroidal field energy and rotational kinetic energy in 
differential motion undergo periodic oscillations on the 



Alfven timescale. The magnitude of the oscillations is in- 
dependent of the initial magnetic field strength; only the 
growth and oscillation timescale depends on the magni- 
tude of the seed field. If viscosity is present, or if some of 
the Alfven waves are allowed to propagate out of the star 
and into an ambient plasma atmosphere, the oscillations 
are damped, rotational energy is dissipated, and the star 
is driven to uniform rotation. 

Recently, Cook, Shapiro, and Stephens ([52], hereafter 
Paper II) generalized Shapiro's calculations for compress- 
ible stars. In their model, the star is idealized as a differ- 
entially rotating, infinite cylinder supported by a poly- 
tropic equation of state. They performed Newtonian 
MHD simulations for differentially rotating stars with 
various polytropic indices and different values of (3. They 
found that when j3 is below the upper (mass-shedding) 
limit for uniform rotation, /3 max , magnetic braking re- 
sults in oscillations of the induced toroidal fields and an- 
gular velocities, and the star pulsates stably. However, 
when (3 exceeds /3 max , their calculations suggest that the 
core contracts significantly, but quasistatically, while the 
outer layers are ejected to large radii to form a wind or 
an ambient disk. 

In this paper, we consider another idealized, but use- 
ful, model to explore the effects of general relativity on 
magnetic braking and viscous damping of differential ro- 
tation in stars. Our star is idealized as an incompressible, 
slowly differentially rotating sphere of uniform density, 
threaded by a seed poloidal magnetic field at time t = 0. 
To simplify the calculations, we restrict our analysis to 
the case in which E Iot <C -Emag <C W, where E I0 t—T 
is the rotational kinetic energy, and E mag is the mag- 
netic energy. The condition E mlig 3> E rot is equivalent 
to the condition that the Alfven timescale tA is much 
shorter than the rotation period P ro t ■ The condition that 
-Emag "C W guarantees that the magnetic field is small 
in comparison to the internal pressure forces and gravi- 
tational field. Wc perform MHD evolution calculations 
in Newtonian gravity, first as a "warm-up" , and then in 
full general relativity, for various initial field configura- 
tions. For an initial magnetic field parallel to the rota- 
tion axis in which we neglect viscosity, the Newtonian 
equations can be solved analytically, but the other cases 
we consider require numerical integration. Adopting ax- 
isymmetry, our MHD equations are two-dimensional, as 
opposed to one-dimensional in Papers I and II. We found 
that this extra degree of freedom changes the description 
of magnetic braking. We found that the angular velocity 
and toroidal magnetic field undergo periodic oscillations 
along the initial poloidal field lines, in the absence of 
viscosity. However, the oscillation frequencies are differ- 
ent along each poloidal field line. The incoherent oscilla- 
tions on different field lines eventually destroys the lam- 
inar flow and creates irregular velocity fields across the 
poloidal field lines. This effect has been studied in New- 
tonian MHD and is sometimes referred to as "phase mix- 
ing" (see [67] and references therein). We demonstrate 
that this phase mixing effect is also present in relativistic 
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MHD. 

In the situation where the Alfven timescale is much 
shorter than the rotation period, which is the limit we 
consider in this paper, the phase mixing may create 
turbulent-like flows in tens of Alfven timescales. We also 
studied the effect of magnetic braking in the presence of 
viscosity. We found that, not surprisingly, the star even- 
tually will be driven to uniform rotation by the combined 
action of the magnetic fields and viscosity. 

Although our analysis in this paper is restricted to 
the early phases of evolution in the slow-rotation, weak- 
magnetic field limit, we are able to solve the full nonlinear 
MHD equations for a highly rclativistic, differentially ro- 
tating star. Our calculations serve as qualitative guides 
and benchmarks for future, more realistic MHD simula- 
tions in full 3+1 general relativity. 

We structure this paper as follows: in Section II, we 
describe our model in detail. We derive and analyze 
the Newtonian MHD equations in Section III. In Sec- 
tion IV, we derive the analogous relativistic MHD equa- 
tions using the 3+1 formalism presented by Baumgarte 
and Shapiro [54] . We then rewrite the MHD equations in 
nondimensional form in Section V, and present our nu- 
merical results in Section VI. Finally, we summarize and 
discuss our results in Section VII. 



II. MODEL 

We consider an incompressible, rotating equilibrium 
star of uniform density (i.e. a polytrope of index n = 0). 
At time t — 0, the star is assumed to rotate slowly and 
differentially with a small axisymmetric poloidal mag- 
netic field. We assume that both the rotational kinetic 
energy and the magnetic energy are small compared to 
the gravitational binding energy. In Newtonian gravity, 
the deviation from the spherical solution is of second or- 
der in the magnetic field strength and/or in the magni- 
tude of angular velocity and is therefore neglected in our 
treatment. In general relativity, the leading order cor- 
rection to the metric comes from the dragging of inertial 
frames, which is of first order. We therefore neglect the 
deformation of the star in Newtonian gravity, but keep 
the frame dragging term in general relativity 

As discussed in Papers I and II, differential rotation 
twists up the frozen-in poloidal magnetic field and gener- 
ates a toroidal field in the star. The toroidal field causes 
a shear stress on the fluid and this stress changes the 
angular velocity profile on the Alfven timescale of the 
poloidal field. Unlike the cases studied in Papers I and II, 
the changes in angular velocity will generate typically a 
meridional current. The interaction of magnetic field and 
the meridional current may develop MHD instabilities 
and result in turbulence [53]. 

In this paper, we consider the case in which the Alfven 
timescale is much shorter than the rotation period. In 
particular, we consider the case in which W ^> E mag 
E mt . The condition E mlig ^ E Iot is equivalent to the 



condition that va 3> £IR, where va is a typical Alfven 
speed, is a typical angular velocity inside the star, and 
R is the radius of the star. It can be shown (see Ap- 
pendix A) that the meridional current can be ignored 
in the early phase of the magnetic braking. In this pa- 
per, we focus on the effect of magnetic braking before 
the meridional current builds up. We assume the system 
remains axisymmetric on the timescale of interest. 



III. NEWTONIAN EQUATIONS 

A. Basic Equations 

We start with the MHD equations for a perfectly con- 
ducting, incompressible, viscous Newtonian fluid (see, 
e.g. [55, 61]) 



dB 

~dt 



= V x (v x B) 



(1) 



B 2 \ 1 



-Vp - pV$ - V — + —(B ■ V)B 

OTT J ATT 



+V • (r/er) + V((V • v) , 
V • v = , 

V 2 $ = AnGp , 

where rj and £ are the coefficients of shear and bulk vis- 
cosity, respectively, and in general are functions of pres- 
sure and temperature. The components of the shear ten- 
sor in Cartesian coordinates are given by 



(2) 
(3) 
(4) 



aij = diVj + djVi - -5 tj V • v 



(5) 



For an incompressible fluid, Eq. (3) allows us to write 



dv 



dB 

— = B ■ Vt) — v ■ VB , 

ot 

1 / B 2 

v-Vv = — Vp-V$-V[ — 

P V 87r 

+V • {va) , 



Aitp 



(6) 



-B VB 



(7) 



where v = r\j p. 

In order to make a direct comparison between New- 
tonian and relativistic MHD equations, we adopt some 
conventions used in relativity For the spherical coor- 
dinate system (r, #,</>), we introduce three orthonormal 
unit basis vectors 

ef = sin 9 cos 4>e x + sin 6 sin <pe y + cos 9e z (8) 

(9) 
(10) 



cos f cos (pe x + cos f sin 0e v — sin Ue z 



ex = -sin0e 2 



cos me v 



where e x , e y and e z are the usual Cartesian unit basis 
vectors. The basis vectors satisfy e~ { ■ e- = 6ij, where i 
and j denote r, 9 and (f>. Any vector V can be expanded 
in these three basis vectors as 



V = V f e f + V e, + V^e, . 



(11) 
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We also define the coordinate basis vectors 

e r = e? , e g = reg, e = rsinfle^ , (12) 

which satisfy • e,- = . Here the spatial metric tensor 
satisfies 

9rr = 1 , 9ee = r 2 , g^ = r 2 sin 2 6 , (13) 

and all the off-diagonal components of gij are zeros. Any 
vector V can be expanded in these coordinate basis vec- 
tors as 

V = V r e r + V e e e + y% , (14) 

hence we have 

V r = V f , V = V § /r, V* = V*/rsm6 . (15) 

We assume the system is axisymmetric and ignore the 
meridional components of velocity (sec Appendix A). 
We therefore set v r — v e — 0, and assume there is no 
toroidal field at time t = 0, B^(0; r,6) = 0. We also 
assume that the star has a small amount of rotation 
= £1(0; r,^) and an initial poloidal field B(0;r,8) — 
B r (0;r,8)e r + B (O;r,9)e e . Equation (6) immediately 
gives dtB r = dtB = 0. Hence the poloidal field, 
which we shall designate -B(O), does not evolve with time. 
Equations (6) and (7) simplify to yield 

d t B^ = B j {Q)d 3 VL , (16) 

d * Q = ~ A T-^-n B °^ d ^ 2 Sin2 9B ^ + (^)™ - 

Anpr 2 sin 8 

where j denotes r and 8, and the usual summation con- 
vention is adopted. In spherical coordinates, this viscos- 
ity term takes the form 

(d t n) vis = la r (^r 4 9 r r») + 1 - d e (v8m 3 6d e Sl) . 
r 4 r 2 sin 9 

(18) 

The toroidal field is generated via Alfvcn waves, which 
cannot propagate in vacuum. Hence the toroidal field 
cannot be carried outside the star. This fact, together 
with Eq. (16), leads to the boundary condition 

B+(t,R,6) = 0=B*(0)d j Sl\ r=R , (19) 
where R is the radius of the star. 



B. Initial Magnetic Field 

In order to solve the MHD equations (16) and (17), 
we need to specify the poloidal field B(0) inside the star. 
In this paper, we consider simple models of the inter- 
nal magnetic field. A systematic way to generate an 
axisymmetric poloidal magnetic field is to assume that 
V x B(0) = inside the star. This corresponds to 
assuming that no internal electromagnetic currents are 



present in the stellar interior. Such a curl-free mag- 
netic field can be generated only by a toroidal current 
on the surface of the star. It follows that we can write 
-B(O) = V$ m , where $ m is a scalar function. The con- 
straint V • B(0) = implies V 2 $ m = 0. The general 
axisymmetric solutions to this equation that are regu- 
lar at the origin can be expanded in terms of Legendre 
polynomials according to 

oo 

<Pm(r,9) =J2 airlpi ( cos ) ' ( 2 °) 

1=0 

where ai are constants. 

The I = mode corresponds to the trivial solution 
B(0) = 0. For I = 1, we have $ m = B rcos9, 
where Bq is a constant. Hence B r (0) = Bo cos 9 and 

B e (0) = —Bosin9, or B(0) = Boe z , which is a uni- 
form field along the rotation axis of the star. It follows 
from Eqs. (16) and (17) that if the angular velocity is 
constant on cylinders, i.e. <9 z il = 0, there will be no mag- 
netic braking in the absence of viscosity. However, it 
can be proven (see, e.g. Ref. [58], Section 4.3) that a 
barotropic star in rotational equilibrium must have such 
a rotation profile. Numerical simulations of core collapse 
supernovae also suggest that the resulting neutron stars 
have angular velocity approximately constant on cylin- 
ders (see e.g., [59]). Hence we also consider the / = 2 
field in which <& m = B r 2 (3 cos 2 8 - 1)/R, where R is the 
radius of the star. The associated magnetic field is 

B f = flo(^) (3 cos 2 9-1) (21) 
B § = -3B sin (9 cos 6> . (22) 

Figure 1 shows the field lines for this I = 2 field. Magnetic 
braking will occur in this case even when d z Sl = 0. As 
will be seen later, the MHD equations for the I = 1 field 
is much simpler than the equations for the I = 2 field. 
For pedagogical purpose, we find it useful to study the 
magnetic braking for the / = 1 field first, where we adopt 
an ad hoc rotation law with d z Q ^ 0. We then study 
magnetic braking for the / = 2 field with a more realistic 
rotation law. 



It is easy to prove that any axisymmetric poloidal field 
can be generated by a vector potential of the form 

A{x) = A\r 1 9)e i . (23) 

For the I = 1 field, we have 

A*(r,6) = ^B rsm9 . (24) 
For the I = 2 field, we have 

A*(r, 9) = B Q sin 9 cos 9 . (25) 
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FIG. 1: Field lines of the I ■■ 
y — (meridional) plane. 
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The vector potential will be useful in Section HID. It is 
also a convenient way to generate axisymmetric poloidal 
fields in relativistic case (see Section IV C). 



C. Conserved Integrals 

The MHD equations (16) and (17) admit two non- 
trivial integrals of motion, one expressing the conser- 
vation of energy and the other conservation of angular 
momentum. 

Multiplying Eq. (17) by pr 2 sin 2 9, integrating over the 
volume of the star and using the boundary condition 
B^(r =R)=0 and the Maxwell equation V • B(0) = 0, 
we find that the angular momentum 



J = j pttr 2 sm 2 9dV 



(26) 



is conserved. 

To derive the energy integral, we first multiply Eq. (16) 
by r 2 sin 2 9B^ /47r and integrate over the volume of the 
star. We obtain, after integration by parts, 



d 
It 



8tt 



dV = —^ j SlBi {Q)dj{r 2 sin 2 OB 4 ') dV 



(27) 

where we have used = rsin9B^, the boundary con- 
dition B^(r = R) = and V • B(Q) = 0. Multiplying 
Eq. (16) by pflr 2 sin 2 9 and integrating over the volume 
of the star, we obtain 



+ J pr 2 sin 2 9 (d t n) vis dV . (28) 
(E mt + £ mag ) = J pflr 2 sin 2 9(d t Q) vis dV , (29) 



Hence we have 
d_ 
~dt 



where 



E„ 



= j ^pfl 2 r 2 sin 2 9 dV 

(B*) 2 



-I 



8tt 



dV . 



(30) 
(31) 



Note that we only include the magnetic energy associated 
with the toroidal field in E mag . One could also include 
the energy associated with the poloidal field and still has 
the energy conservation since the poloidal field does not 
change with time. If we define 



where 



we have 



E v 



E vis (t) - / E vis (t')dt' , 
Jo 

, = J pflr 2 sm 2 9(d t n) vis dV 



+ E V 



. 



Hence the total energy E — E mag + E mag + E, 
served. 

To compute E v i s , we use Eq. (18) for (d t fl) 
integration by parts, we obtain 



(32) 

(33) 

(34) 
vis is con- 
After 



E v 



I 
I 



rjr sin 9 



(d r nf + 



j]r 2 sin 2 6»(VO) 2 dV . 



dV , (35) 
(36) 



i j \ P ^r 2 ^ 2 9dV = ±J Bi(0) dj (r 2 sin 2 0B*)dV ^ = Bj{Q)d . n 



We note that, consistent with the nonrelativistic MHD 
approximation, the electric field energy E 2 /8ir is not in- 
cluded in Eq. (34) and the angular momentum of the 
electromagnetic field is not included in Eq. (26) [56]. 

The motivation for monitoring the conservation equa- 
tions during the evolution is twofold: physically, evalu- 
ating the individual terms enables us to track how the 
initial rotational energy and angular momentum in the 
fluid are transformed and/or dissipated; computation- 
ally, monitoring how well the conservation equations are 
satisfied provides a check on the numerical integration 
scheme. 



D. Build up of Small-Scale Structure: Phase 
Mixing 

In the absence of viscosity, the MHD equations (16) 
and (17) become 



(37) 
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Anpr 2 sin 



-B j (0)dj(r 2 sin 2 OB*) . (38) 



These equations have a peculiar property that will in- 
duce the growth of small-scale fluctuations in the fluid's 
angular motion and toroidal field. To study this, we first 
notice that the 2+1 MHD equations can be reduced to a 
1+1 problem if we can find a coordinate transformation 
u = u(r, 9), v =v(r, 9) such that 



B^O) 



d 



k(u, v) 



d_ 

du 



(39) 



If such coordinates exist, the MHD equations (37) 
and (38) reduce to the following 1+1 form: 



d t B* 



k(u, v)d u fl , 
1 

4irpr 2 sin 2 9 



(40) 

K(u,v)d u {r 2 sin 2 6 B+) . (41) 



Here r and 9 are treated as implicit functions of u and v. 
Applying Eq. (39) to u and v, we obtain the conditions 
that must be satisfied by u and v: 



= o . 



(42) 
(43) 



Such coordinates indeed exist and are easy to find. 
Equation (43) requires that v is constant along a given 
poloidal field line. Since the field lines do not cross and 
cover the volume insider the star, v is a valid coordinate, 
labeling each field line. One way to find v is by way of 
the vector potential. The covariant ^-component of the 
vector potential is defined as 



rsm 



OA* . 



It follows from B(0) = Vx {A*e A ) that 

B r (0) 
B e {0) 



r 2 sin 2 9 
r 2 sin 2 9 



(44) 

(45) 
(46) 



It follows that B 1 * (fydjAtj, = 0. Hence we can choose 
v = f(A < / > ), where / is any function with a continuous 
first derivative. The coordinate u is arbitrary as long as 
it is independent of v. The function n is then determined 
by Eq. (42). 

It is obvious that for the / = 1 field (B(0) = Boe z ), 
the preferred coordinates are cylindrical coordinates (u = 
z = rcos9 and v = zu = rsin9), for which n — B n . The 
MHD equations become 



d t B* - B Q d z n , 
d t n = ^d z B* , 



(47) 
(48) 



which can be solved analytically (see the next subsec- 
tion). 

The situation for the I = 2 field is more complicated. 
The covariant ^-component of the vector potential that 
generates the field is A^ = B (r 3 /R) sin 2 OcosO [see 
Eq. (25)]. We choose 



u = -r 2 (3cos 2 6»- 1) , 
v = r 3 sin 2 9 cos 9 . 



(49) 
(50) 



The variable u is chosen so that Vm cx B(0), hence u cx 
$ m , and Vti • Vd = 0. Equation (42) gives 



2B r 2 
R 



(3 cos 2 9 + 1) . 



(51) 



The MHD equations become 



d t B* 



2B r 2 

R 
2Bo_ 
AirpR 



(3cos 2 9+l)d u n (52) 
[r^Scos^ + l^B'^ + .B^ .(53) 



Here r and 9 are regarded as implicit functions of u and 
v. 

The set of equations (40) and (41) is effectively a 1+1 
system, with v serving as an effective "mode" number. 
There is no coupling between different values of v. The 
toroidal field B* and angular velocity £1 oscillate inde- 
pendently along each poloidal field line. The oscillation 
frequencies are determined by the boundary condition 
B'f'ir = R) = and are in general different on differ- 
ent field lines. As a result, the toroidal field and angular 
velocity on different poloidal field lines will lose coher- 
ence. Small-scale structure will gradually build up. This 
is commonly referred to phase mixing [67, 68]. In the 
next subsection, we demonstrate this phenomenon quan- 
titatively via an analytic solution. 



E. Analytic Solution for the 1 = 1 field 

Combining Eqs. (47) and (48), we obtain 
,2 U 4> _ v 2 A d 2 B' t ' , 



dfBf 



(54) 



where va — B n / y/Anp is the Alfven velocity associated 
with the poloidal field. This is a simple wave equation; 
a similar equation arises for f2. We are interested in the 
solution in which f2 is symmetric under z — > — z and 
B^ = at t = 0. The general solution that satisfies the 
boundary condition B^(r = R) = is 



oo 

B*(t;w,z) = B (^^C n {w)wi[k n (w)z]i 



(55) 
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fl(t;zu,z) — £7(0; zu, z) + C n (ru) cos[k n (m)z} x 

{1 - cos[w„(ro)i]} , 



where 

k n (m) 



rnr 



VR 2 - tu 2 



(56) 

, w n (ro) = «AA:„(tn) . (57) 



The functions C n {w) are determined by the initial dif- 
ferential angular velocity distribution Q(0:zu, z): 



2R 



uj n \jR 2 — uj 2 Jo 



•JW- 



[<9 z fi(0; w, z)\ sin k n zdz . 

(58) 

We see clearly from Eq. (57) that the characteristic 
frequency on each field line is different. Incoherence and 
small-scale structure will build up after a certain time. 
To be precise, for a given w and a small length scale AL, 
the eigenfrequency to n (w) changes on this length scale by 
an amount 



Aw n 



(R 2 - vo 2 fl- 



-AL . 



(59) 



The field and angular velocity will lose coherence 
over AL on a timescale 



2tt 



2R 2 



.... , . ^ 3/2 

*A^„ n-AI.V R 2 



tA , 



(60) 



where the Alfven time tA = R/va- For w = R/2, AL = 
R/10 and n = 1, we have £ co h ~ 26^. The phase mixing 
will eventually lead to chaotic-like motion over any finite 
radial scale. 

The solution for the 1 = 2 case must be obtained nu- 
merically. We postpone a discussion of this case, and 
the effects of viscous damping, for the fully relativistic 
treatment. 



IV. GENERAL RELATIVISTIC MHD 
A. 3+1 Decomposition 

In this subsection, we briefly summarize the gen- 
eral relativistic MHD formulation presented in Ref. [54]. 
Hereafter, we adopt geometrized units and set c = G = 1. 

We start with the standard 3+1 decomposition of the 
metric: 

ds 2 = -a 2 dt 2 + ll0 {dx l + f3 l dt)(dx ] + (3 j dt) . (61) 
The spatial metric j ab is related to the full metric g a b by 
Jab = 9ab + n a n b , (62) 



where n a is the unit normal vector n a = —aSI a t to the 
spatial slices. Next we decompose the Faraday tensor as 



F ab = n a E b _ n b E a + £ o6c r 



(63) 



so that E a and B a are the electric and magnetic fields 
observed by a normal observer n a . Both fields are purely 
spatial, i.e. 

E a n a = = B a n a . (64) 
The three-dimensional Levi-Civita tensor is defined by 

(65) 



^a6c 

abed 



e abcd n d 



1 



[abed] , 



(66) 



where [abed] is the antisymmetrization symbol. We 
also decompose the current four-vector J a as 



J a = n a Pe + J a 



(67) 



where p e and J a are the charge and current density as 
observed by a normal observer n a . Note that J a is purely 
spatial. 

With these definitions, Maxwell's equations VbF ab = 
4irJ a and V[ a F bc ] = take the following 3+1 form: 

D i E i = 47r Pe , (68) 
dtE i = t Hk D .( aBk ) _ 47ra j* + aKE i + £ E i , (69) 

D i B i = , (70) 
8 t B l = -e ijk Dj(aE k ) + aKB 1 + £ p B l , (71) 

where K is the trace of the extrinsic curvature, Dj is the 
covariant derivative compatible with 7^ , and £ denotes 
a Lie derivative. 

It follows from V[ a F bc } = that F a b can be expressed 
in terms of a vector potential 



Fab = d a A b - d b A a . 
We decompose A a according to 

A a = ^cTla + A a , 



(72) 



(73) 



where the three- vector potential A a as observed by a nor- 
mal observer n a is purely spatial (A a n a = 0). The rela- 
tionship between $ e , A^, E % and B^ are given by 

d t Ai = -aEi - di{a$ e ) + £fsAi , (74) 
B* = e ijk A ktj . (75) 

In the ideal MHD limit, the fluid is assumed to have 
perfect conductivity, which is equivalent to the condition 
that the electric field vanish in the fluid's rest frame. The 
relation between the electric and magnetic field is given 
by 



aE t 



(vi+f3i)B k , 



(76) 



8 



where v J = vP /it*. In the Newtonian limit, the above 
equation reduces to the familiar form E = —v x B. 

In the ideal MHD limit, the Faraday's law (71) be- 
comes 



d t B l = d j {v i B j -v'B^ , 



(77) 



where B l = ^/jB l . Here 7 is the determinant of the spa- 
tial metric 7^ . Equation (77) is the relativistic version 
ofEq. (6). 

To derive the relativistic Navier-Stokes equation in the 
ideal MHD limit, we separate the stress-energy tensor 
into a fluid part and an electromagnetic part: 



nab 



rpab 1 rpab 
1 fluid ' 1 cm 



where 



rpab 
1 fluid 

nab 



rpab _j_ rpab 
1 p ~r 1 vis 



T p a0 = phu a u° + Pg 



^ab 



rpab 
em 



_L I pacpb _ I„afcp pcd 

4tt *v c 4 5 cd 



(78) 

(79) 
(80) 

(81) 



Here /i = 1 + e + P/pis the specific enthalpy, pe is the 
internal energy density, and Tjf b is the stress-energy ten- 
sor for a perfect fluid. An incompressible fluid may be 
regarded as a gamma-law equation of state P = (T — l)pe 
in the limit r — > 00. Hence we have e = for an incom- 
pressible fluid. The viscous stress-energy tensor T" b s is 
given by [see, e.g. [60], Eq. (22.16a)] 



rpab 
vis 



2r]a ab - (P ab V c u c . 



(82) 



The projection tensor P ab and the shear tensor a ab are 
defined as 



pab = g ab + u a u b 

<j ab = i(V c u a P cb + \7 c u b P ca ) 



(83) 

-^P ab X7 c u c . (84) 



For an incompressible fluid, V c u c = 0. The relativis- 
tic Navier-Stokes equation is computed from the spatial 
components of the equation V ' a T ab — 0. The result is 

dt(VlS t ) + d^Sivi) = -a^y (&iP + 

+a^F la J a + 2d b (a^ W b ) + Wbc d l9 bc , 

where 

S a = phau l u a . 
Lorentz's force law implies 



47^7° = V b F ab 
Hence we have 



1 



ab\ 



a^F ta J a = — F la d b {a^jF ab ) 



(85) 
(86) 

(87) 
(88) 



Equations (85) and (88) comprise the relativistic version 
of Eq. (7). 



B. Basic Equations 

The metric for a uniform density, spherical star is given 
by (see, e.g. [60], Box 23.2) 



ds 2 = 
where 

a(r) 
and 



-a 2 (r)dt 2 



dr 2 
A2(r) 



r 2 {d9 2 + sin 2 



! ), (89) 



A(r) 



i-¥ -h^- 2 -^ (r<R) 



2M 



(r > R) 



1 - 



2Mr 2 
R 3 



The interior pressure is 

A(r) - a R 



(r < R) 
(r > R) 



P(r) 



3a R -\{ry 



where 



2M 

a R = a (r = R) = \jl- — 



(90) 



(91) 



(92) 



(93) 



In the presence of slow rotation, small magnetic field 
and viscosity, the diagonal components of the background 
metric has a correction of order fi 2 , |£?| 2 and r/|Vf2|, 
which will be neglected since they are of higher nonlinear 
order. Rotation induces the off-diagonal component g t< p, 
which corresponds to the dragging of inertia frames. We 
write 



5t0 = P<j> = -r 2 sin 2 9 tj{r, 9) 



(94) 



where w (r, 9) is the angular frequency of the zero angu- 
lar momentum observer (ZAMO) observed by an inertial 
observer at infinity. In order of magnitude, lu/Q ~ M/R. 
We adopt the Cowling approximation, which assumes 
that the background metric is fixed, i.e. <9t<? M „ = 0. We 
can justify this approximation by the following argument. 
First, we can use the gauge degrees of freedom to freeze 
the lapse and shift: dtot — dtfti — 0. Second, the evolu- 
tion of the spatial metric 7^ is given by the ADM equa- 
tions [69] 

d t jij = -2aK l3 + Drfj + Djf3i , (95) 
d t Kij = -D l D J a + a{R lJ -2K lk K k J +KK lJ ) 

1 



-8ira 



Sij 



-,lij(S-Ps) 



+K lk D 3 (3 k + K k3 Di(3 k 



+ l3 k D k K l] 

(96) 

where K~ij is the extrinsic curvature, Rij is the three- 
dimensional Ricci tensor, K = 

p s = n^n v T^ v , 

Sij — "iik^ijmF 



K^j, and 



S = S j , 



(97) 
(98) 
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where T^" is the total stress-energy tensor. We assume 
that the star is in hydrostatic quasi-cquilibrium at t — 
0. The time derivatives dtjij and dtKij come from the 
changes in the matter and B- field sources term Sp s , 5Sij 
and SS, which are all of the order fi 2 , \B\ 2 and ^|VQ|. 
This means that the deviation of jij from its initial value 
will remain higher order. Hence our gauge conditions 
does not introduce a first order correction to the spatial 
metric in later times, and the Cowling approximation is 
valid. 

As in the Newtonian case, we assume that the sys- 
tem is axisymmetric and v r (t;r,9) — v e (t:r 7 9) = 0, 
^O;r,0) = n(t;r,6), 5^(0; r,6) = 0. It follows from 
Eq. (77) that d t B r = d t B = 0. Hence the poloidal field 
remains unchanged, just like the Newtonian case. Denote 
the initial magnetic field 

B(0; r, 6) = B r (0; r, 9)e r + B e (0; r, 9)e g , (99) 

where e r — d/dr and eg = d/d9. Taking the <fi- 
componcnt of Eq. (77), we have 



d t B* = J-djiSlyftB*) 

= n-^dji^B^ + B^jn 

= B 3 {<d)d 3 n , 



(100) 



where we have used Eq. (70) to obtain the last equality. 
The 0-components of Eqs. (85) and (88) simplify to 



d t S^ 



1 



F <pb d c (a^fF bc ) + 2d b (a y ^ V al) . (101) where 



471-/7 

It follows from = fi and u a u a = —1 that 



u l = - + o(n 2 

a 



= - + o(n 3 ) 

a 



(102) 



r 2 sin 2 i 



a 



■(n-uj) + 0(n 3 ) . (103) 



Using Eq. (86), we obtain 



P + P r 2 sin 2 0(n-u) + 0(n 3 ) . (104) 



2<) l Ao x - W b <f > ) = d r {r]Xr 4 sin 3 6d r n) + 8 e ( r f ^ 6 d e n 



+o( v n 2 ) . 

Gathering all the terms, we obtain 



(106) 



1 + 



7 «£*(0)fl'(0) 



4tt(p + P) 
+ 



d t n-- 



(p + P)r 4 sin 3 9 



OiB? {0)dj{ar 2 sin 2 
4tt(p + P)r 2 sin 2 9 

[<9 r (r/r 4 Asin 3 9d r U) 



r 2 sin 3 i 



A 



(107) 



The term on the left-hand side, j ij B i (0)B^ {0)/An{p ■ 



P) 



E ma .g/M, is assumed to be small. Hence it 



can be neglected. We also note that to the leading order 
in our adopted gauge, the frame dragging frequency w 
does not enter into the evolution equations. It follows 
from Eq. (92) that 



p + P(r) = p 



an 
a(r) 



(108) 



Combining the results, we obtain the relativistic version 
of Eqs. (16) and (17): 



dtB 1 * 



H d = B J {0)d,n 

no a " 1 
d t il = 



^-B j (0)di(ar 2 sin 2 9B" 
a R Airpr 2 sm 2 9 ' JV 

+(d t n) vis , 



(109) 



(110) 



1 



Ar 2 sin 3 9 



de(vsin 3 9d n) 
(111) 



a 



Here we again set v = rj/p. The condition that the 
toroidal field cannot be carried outside the star gives the 
boundary conditions B^ = at r — R, which also implies 
[from Eq. (109)] B^(0)djfl = at r = R. In the New- 
tonian limit, with M/R < 1, Eqs. (109)-(111) reduce to 
the Newtonian MHD equations (16)-(18). 



C. Initial Magnetic Field 



Straightforward calculations from Eqs. (63) and (76) us- 
ing the metric (89) yield, to the leading order, 



F^ b d c (a^F bc ) 



r 2 sin ( 



B j (0)dj(ar 2 sin 2 0B 4 



r 4 sin 3 9 

.L±_l lijB * { o )B i(o)d t n , (105) 



where we have imposed our gauge condition d t Pi = 
by dropping a term involving dtOJ- The viscosity term is 
[from Eq. (84)] 



The two sets of poloidal field functions we discussed in 
Section IIIB cannot be used in the relativistic case be- 
cause they do not satisfy the relativistic Maxwell equa- 
tion DiB l ({)) = 0. We want to choose two sets of initial 
field which will reduce to the ones we discussed in Sec- 
tion IIIB in the Newtonian limit. It is easy to prove that 
any axisymmetric poloidal field can be generated by a 
vector potential of the form Ai = A^Sf . Hence we can 
specify A^, and then compute £? J (0) from Eq. (75). This 
guarantees that the constraint equation DiB l (0) = is 
satisfied. 
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The simplest way to generalize the Newtonian initial 
fields is to use the same as in the Newtonian cases 
[Eqs. (24) and (25)]. Hence the "I = 1" field is generated 
by 



M = \b t 2 sin 2 9 . 



(112) 



The corresponding poloidal field is given by Eq. (75), 
yielding 



B, 



B r = XB cos9 , B° = -A— sin 9 . (113) 

r 

We define cylindrical coordinates vo = r sin#, z = rcosO. 
The cylindrical components of the poloidal field are 



B™ — , B z = XBn 



(114) 



We see that the poloidal field is still in the z direction. 
However, its amplitude 



B(0)\ = JjijB^BiiO) 



2Mw 2 



B \ 1 



R 3 



(115) 



decreases with increasing vo. 

The generalized u l = 2" field is generated by 

r 3 

= B — sin 2 9 cos 9 . 
R 

The poloidal field has components 



(116) 



2 ° ^ ° e - sin 9 cos 9 . 



B r = \B -(3cos z 9-l) , B' 
R 



R 



(117) 

ft follows that the shape of the field line of this I = 2 
field is the same as the Newtonian I = 2 field shown in 
Fig. 1, even though the magnitude is different. 



D. Conserved Integrals 

As in the Newtonian case, the relativistic MHD equa- 
tions (109) and (110) admit two non-trivial integrals of 
motion, expressing the conservation of energy and an- 
gular momentum. The derivation is very similar to the 
Newtonian case, so we just state the result here. 

The conserved energy contains a sum of three terms: 



where 

-Erot 



B -EVot ~i" B m ag ~t" E v [ s , 



P ^ B 2 2 • 2/1 

\l r sin 



(118) 



-drd6d(j), (119) 



E, 



mag 



{B^ f r 2 sint 
8?r A 



drdOdcj) 



E v 



f E vis (t')dt' 
Jo 



(120) 
(121) 



and where 

E vis = J r]r 2 sm 2 6 



x 2 (d r n) 2 + 



r 2 sin t 



= J rjr 2 sin 2 6g ij ViQVjCl 



r 2 sin 



A 

-drdOdcj) . 



drdOdcj) (122) 
(123) 



The conserved angular momentum is 



J ol 



P 



-{Ur 2 sin 2 9) 



r sm ( 



drd6d<j> . (124) 



E. Special Coordinates and Phase Mixing 

As in the Newtonian case, in the absence of viscosity, 
Eqs. (109) and (110) can be brought into 1+1 form by 
the coordinate transformation u = u(r,6), v = v(r,9) 
that satisfies 



B'(«)^j = «(«,«), 
*(„>*- 0. 



(125) 
(126) 



In these new coordinates, Eqs (109) and (110) become 
(without viscosity) 



dtB* = nd u Q. 



(127) 

d u (ar 2 sin 2 OB*) . (128) 



an Anpr 2 sin 2 9 
It follows from Eq. (75) that when Ai = A^Si 

B r = -±-d e A^ , B e = -^-d r A^ . 



(129) 



Hence we have B t diA c f > — 0. We can choose v — f{A$) as 
in the Newtonian case, where / is an arbitrary function 
with continuous first derivative. Since we use the same 
A^ as the Newtonian case in constructing the initial mag- 
netic field, we can use the same u and v to simplify the 
MHD equations. 

In particular, we use cylindrical coordinates defined by 
zn = rsin# and z = rcos9 for the I = 1 field. In these 
coordinates, the MHD equations become 

dtB 4, = XB d z Q (130) 

dt n = ^ d z (aB*) • (131) 

a R \^pj 



For the / = 2 field, we choose 



u = -r 2 (3cos 2 6»- 1) , 
v = r 3 sin 2 9 cos 9 . 



(132) 
(133) 



It follows from Eq. (125) that 
K = B j (0)d jU -- 



2AS ° r2 (3cos 2 0+l) . (134) 



R 
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The MHD equations become 
^ = _2A|o^ (3cos2( 



R 

—)^%[r 2 (3cos 2 6 + l)d u (aB*) 
a R J AirpR' 



(135) 



+aB 4 '} 



(136) 



Here r and 9 are regarded as implicit functions of u and 
v. 

As discussed in Section HID, Eqs. (127) and (128) im- 
ply that there is no coupling of and f2 between dif- 
ferent poloidal field lines. The toroidal field and angular 
velocity will oscillate in general at different characteristic 
frequencies along different field lines. The angular veloc- 
ity flow and toroidal field pattern will become irregular 
as a result. We will demonstrate this phase mixing effect 
in Section VI when we present our numerical results. 



V. NONDIMENSIONAL FORMULATION 

Before numerically integrating the MHD equations, it 
is convenient to introduce the following nondimensional 
variables: 



r = 



B = 



B j (0) 



v = — 



E 
J 




(137) 
(138) 

(139) 

(140) 
(141) 
(142) 
(143) 
(144) 

(145) 



Here fi c is an arbitrary constant with a magnitude char- 
acteristic of the initial angular velocity. In terms of these 
new variables, Eqs. (109) and (110) become 



M{m c ) 2 
J 



MR 2 n c ' 

where the Alfven time tA is defined as 

R 

tA = 



d { B = aB 3 {0)d-tt 
a 2 &{<)) 



(146) 



i" - — ^^^(f 2 sin 2 ^) + (9^) vis , (147) 
olr r z sin J 



where <9j = Rdj and 



(0 f fi), 



a#f 4 sin 3 9 



sin' 



dr{vr 4 sin 3 OOfST) 

3/3 



-dt>n 



(148) 



Since the background pressure distribution is spherical, 
it is reasonable to assume that v = v(r), hence 



OiR 



df{v\) + -v\ 



+ -^{d 2 e Cl + 3 cot Odgfl) 



}' 



(149) 



The Navier-Stokes equation and the regularity condition 
on the surface of the star require that the dynamic vis- 
cosity r] in general has to decrease to zero continuously 
as the surface of the star is approached. We adopt the 
simplest viscosity model in our numerical computations: 
we assume that v = r\j p is essentially constant in the 
interior, but rapidly drops to zero when approaching the 
surface of the star at r = R. 

While the variable B is convenient for numerical evo- 
lution of the MHD equations, we also found it convenient 
to introduce another non-dimensional variable 



B 



f sin 9 




(150) 



which is a measure of the strength of the toroidal mag- 
netic field. Note that in this non-dimensional formu- 
lation, we do not need to specify the ratio va/Q c R, 
which is assumed to be much greater than one so that 
Em ag ^> E TOt . We also do not need to specify the char- 
acteristic amplitude of the seed poloidal magnetic field 
B , which has to be small so that E mag <C W or cquiv- 
alently, va -C 1. These quantities are absorbed in our 
nondimensional variables. 

In the next two subsections, we will further simplify the 
general equations above by inserting the specific forms of 
the two sets of initial magnetic field. 



A. Equations for the 1 = 1 Field 

As discussed in Sections III D and IV E, the cylindrical 
coordinates w = f sin 9 and z = f cos 9 are the preferred 
coordinates in the absence of viscosity. To impose the 
boundary condition at f = 1, it is convenient to intro- 
duce a new variable u = z/yl — m 2 in place of z. The 
boundary of the star is located at it = 1. We are inter- 
ested in a solution in which f2 is symmetric under the 
reflection z — > —z. Hence the computation domain is 
u e [0, 1] and w G [0, 1]. The MHD equations, in terms 
of u and w, take the form 



diB = 



vT 



(151) 



12 



dp 



a 2 X 



- T duB + (dtti), 



(152) 



where (<9 t -f2) v is is obtained from Eq. (149) and Eqs. (C7)- 
(CIO). 

The nondimcnsional energy and angular momentum 
are given by 

J = 3or / rfwVl-w 2 / du— (153) 
Jo Jo a 2 A 



£ = E Iot (t)+E, 



mag(f) + / 

Jo 



^vis(*') , (154) 



where 



E mt = I [ dmw^l-w 2 e\°\w) 
2 Jo 



Er, 



E v 



dww z y/l — w 2 e™ 3,8 (w) 



\ [ dm [ 
Jo Jo 



1 

dz 

A 



(155) 
(156) 
(157) 



The functions e[ ot (w) and e™ ag (ro) are given by 



er s (^) = f 

Jo 



du- 



a 2n2 

e 2 A 

B 2 



du^r-O, , 

a^A 



A 



(158) 
(159) 



In the absence of viscosity, there is no coupling between 
different values of w. Hence the reduced energy and an- 
gular momentum functions 



ji(w) = f du-^- 

Jo Oi * 



(160) 
(161) 



are also conserved for each value of w. This can also be 
proved directly from Eqs. (151) and (152). 



B. Equations for the I = 2 Field 

In the absence of viscosity, the preferred coordinates 
are u = — r 2 (3cos 2 6* — 1) and v = r 3 sin 2 cos 0. In 
numerical calculations, it is more convenient to introduce 
the following nondimcnsional variables 



ui(v) 



u 2 (v) - Ui(v) 



3V3 



r 3 sin 2 6 cos ( 



(162) 
(163) 



where ui(v) and u 2 (v) are the values of u at which a 
given v=constant line intercepts the sphere r = R, the 
surface of the star, with u 2 > U\ (see Fig. 2). 




s '% t \u 2 (v) 







FIG. 2: The points u\(v) and U2(v) 



The MHD equations, in terms of u and v, take the 
form 



di B = _^Ar 2 (3cos 2 ^l)^ 



u 2 (v) — u\(v) 



dp 



2a 2 \ 



! (3cos 2 



u 2 {v) - ui(v) 



duB + B 



(164) 



(165) 



Note that the relativistic factors a and A are functions 
of r. The variables r and 9 are regarded as implicit func- 
tions of u and v. The transformation between these two 
sets of variables are derived in Appendix B. we only con- 
sider the solution with equatorial symmetry. The compu- 
tation domain is u € [0, 1] and v e [0, 1]. The boundary 
of the star is located at u — and u = 1. As in the New- 
tonian case, lines of constant v coincide with poloidal 
field lines. 

It is not easy to handle the viscosity term (dp) v - ls nu- 
merically in this (u, v) coordinate system because of the 
coordinate singularity at v = and = 1. Besides, these 
particular coordinates lose their advantage in the pres- 
ence of viscosity. Hence we use the cylindrical-like co- 
ordinates introduced in Section VA when we study the 
effect of viscosity. 

In the absence of viscosity, the energy E and angular 
momentum J are given by 



where 



E 
J 

e 2 (v) 



— r. 



u 2 {v) - ui(v) „ 



2V3 7 di R 2 m (166) 

T3 J o dv - J2 (v), (167) 



;rot 



e r 2 ot (0) = / du— 
Jo a ' 



an sin 



A(3cos 2 6> + 1) 



Q 2 



(168) 
(169) 
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TABLE I: Summary of the cases we have studied. 



Case 


Gravitation 


M/R 


Poloidal field 


Viscosity 


la 


Newtonian 


< 1 


I = 1 





lb 


GR 


0.3 


I = 1 





Ic 


GR 


0.44 


/ = 1 





Ha 


Newtonian 


< 1 


/ = 2 





lib 


GR 


0.3 


/ = 2 





lie 


GR 


0.44 


/ = 2 





III 


GR 


0.3 


I = 2 


0.002 



-mag 



(*) 



du 



sin 



hi 



v) = du— 
Jo " 



a\(3cos 2 9 + 1) 
ft sin 2 6 



2 A(3cos 2 6» + l) ' 



B 2 , (170) 



(171) 



It follows from Eqs. (164) and (165) that in the absence 
of viscosity, the functions e 2 (v) and j2(u) are conserved 
for each v. 
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VI. NUMERICAL RESULTS 

We numerically integrate the two sets of differential 
equations (151), (152) and (164), (165) by the iterated 
Crank-Nicholson method [57]. A more detailed descrip- 
tion of our finite differencing scheme is given in Ap- 
pendix C. In this section, we report our numerical re- 
sults. Tabic I summaries the cases we have studied. 



FIG. 3: The 6 = 45° line (line A) and vj/R = 0.5 line (line 
B). The circular arc is the surface of the star r — R. 



in Eq. (58). Hence in the Newtonian limit (M/R < 1), 
the toroidal field and angular velocity oscillate with the 
fundamental cigcnfrcqucncy on each cylindrical surface, 
given by Eq. (57) for n = 1. The analytic solution in this 
limit is 



A. Without Viscosity (Cases la— lie) 

In the absence of viscosity, we use the special coordi- 
nates discussed in Sections V A and V B to integrate the 
nondimensional MHD equations. 



1. 1 = 1 initial field (Cases Ia-Ic) 

As discussed in the Newtonian case, there is no mag- 
netic braking if the initial angular velocity distribution is 
constant on cylinders for the I = 1 initial field. The same 
result applies for the relativistic case, as is indicated from 
Eq. (130). To study magnetic braking, we use the follow- 
ing ad hoc initial differential angular velocity profile [70] 



O(0; w, z) 



I 



l-w 2 



1 



VT^W 2 



(172) 



which is not constant on cylinders. This corresponds to 
setting the function 



C n (w) 



Iva V TP) 



n = 1 
n > 1 



(173) 



fo(t;zb,z) = 



1-w 2 



l + w 2 2 
cos 



COS 



7TZ 



VI- w 2 , 



B(i;w,z) = 



\-w 2 



sin 



7TZ 



Vi- w 2 



7rf 



,Vl - W 2 
It follows from Eqs. (158)-(161) that 



1 il-vj 2 ) 2 
+ ^ — — x 



-mag / ? - \ 



(l + w 2 ) 



cos 



irt 



Vl - W' 



,2\2 



- — sin 2 



Tit 



vr^ 2 



e\(w) 
hi™) 



(\+VJ 2 ) 2 
1 



l+VJ 2 ' 



(174) 



(175) 



(176) 



, (177) 



1 (178) 



(179) 
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FIG. 4: Snapshots of B (solid lines) and Q (dashed lines) 
along the 9 = 45° line (line A of Fig. 3) for a Newtonian star 
(M/R <C 1). The initial magnetic field is the I — 1 field given 
by Eq. (113). As expected, small-scale structures build up at 
late times. 
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FIG. 5: Snapshots of B (lower curves) and Q, (upper curves) 
along the w/R — 0.5 line (line B in Fig. 3) for a Newtonian 
star (M/R <C 1). The surface of the star is located at z = 
zr — \/R 2 — vj 2 — R/\/2. Each curve represents the profile 
at a particular time. Both B and tl oscillate sinusoidally with 
a period T — \fi = 1.73. Curves are plotted at time intervals 
At = 0.2 from t = to t = 8.6 (about 5 oscillation periods). 



We plot the snapshots of B and Cl as a function of 
time along two arbitrary lines as shown in Fig. 3. Fig- 
ure 4 shows the snapshots of B and fi as a function of 
time along the 9 = 45° line based on the analytic solution 
given by Eqs. (174) and (175). We see that small-scale 
structures develop at late times, as predicted from our 
previous analytic study. Figure 5 shows the snapshots 
along the w/R = 0.5 line. We see simple sinusoidal oscil- 
lations with a period T = -\/3 = 1.73. No small structure 
is observed since we are looking in the direction along a 
particular poloidal field line. 

Since j\ (w) is conserved, we can define a mean angular 
velocity for each value of w by 



where 



h{w) 



r 1 du 

Jo 



In our Newtonian example, we have 



1 



l + m 2 



(180) 



(181) 



(182) 



Equation (174) indicates that for each w, the angular 
velocity Cl oscillates about Ct mean with a period T = 
2-s/l - w 2 - The first term of Eq. (176) is the reduced 
rotational kinetic energy associated with f2 mca n, which is 
independent of time. The kinetic energy associated with 
the differential rotation transfers back and froth with the 
energy associated with the toroidal magnetic field with a 
period T/2 (see Fig. 6). 



0.05- 



-0.05 




FIG. 6: Energy evolution for a differentially rotating New- 
tonian star with an I = 1 poloidal magnetic field. The solid 
line shows e\ ot (i) - e\ ot (0) , the dotted line shows e™ s (i). The 
dashed line is the sum of the two lines, which is ei — e\ ot (0) = 
All the reduced energies are evaluated at w = 0.5. 



There is no analytic solution to the MHD equations 
when the star is relativistic. We solve the MHD equa- 
tions (151) and (164) numerically in this case. The ap- 
pearance of small-scale structures does not cause any 
problem in our numerical code, because we have chosen 
suitable coordinates to avoid taking derivatives perpen- 
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FIG. 7: Same as Fig. 4 but for a relativistic star with M/R — 
0.3. The behavior is qualitatively the same as the Newtonian 
case, but the amplitude of oscillation is larger. 



dicular to the poloidal field lines. To check our numerical 
integrations, we perform second-order convergence tests 
and, in addition, monitor the conserved integrals e\{w) 
and ji(in) defined in Eqs. (160) and (161). We find that 
the fractional change of i\ to be less than 10 -4 from time 
t = to t = 100 for various values of the compaction ra- 
tio M/R. The fractional change of j\ is less than 10~ 5 
over the same period of time. 

Figure 7 shows the snapshots of B and £1 for a rel- 
ativistic star with M/R = 0.3 along the 9 — 45° line. 
We see that the behavior is qualitatively the same as the 
Newtonian case, but the amplitude of the toroidal field 
is larger by a factor of 2. Figure 8 shows the snapshots 
along the vj / R — 0.5 coordinate line. The initial angular 
velocity distribution [Eq. (172)] no longer corresponds to 
the fundamental mode of the relativistic equations, but 
it is still very close. The oscillation frequency is found to 
be T = 4.3, which is longer than the Newtonian period 
by a factor of 2.5. The corresponding reduced energy is 
similar to Fig. 6, but the amplitudes are larger and the 
period is longer by a factor of 2.5 (Fig. 9). This can be 
understood by rewriting Eqs. (130) and (131) as 



OtR 



djB + 



d £ (a 2 X) 
a 2 X 



B 



(183) 



If the background metric varies slowly along the field 
lines so that the term dz(a 2 X)/a 2 \ <C 1, the equation 
simplifies to 



Q. 0.8 



B 
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FIG. 8: Snapshots of B (lower curves) and Q, (upper curves) 
along the zu/R — 0.5 line for a relativistic star with M/R — 
0.3. As in the Newtonian case, both B and Cl oscillate period- 
ically. The oscillation period is T = 4.3, which is longer than 
the Newtonian case because of the relativistic time dilation 
(see the text). Curves are plotted at time intervals At = 0.4 
from t = to t — 21.2 (about 5 oscillation periods). 




FIG. 9: Same as Fig. 6 but for a relativistic star with M/R = 
0.3. 



where 



a 3 \ 2 

OLR 



1/2 



t 



(185) 



dlB w dlB 



(184) 



Note that Eq. (184) does not contain explicitly the rela- 
tivistic factors a and A, which are absorbed in the new 
time variable f and B [see Eq. (150)]. Hence in this ap- 
proximation, the relativistic corrections are that (a) time 
is dilated according to Eq. (185), and (b) the toroidal field 
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FIG. 10: Same as Fig. 4 but for a relativistic star with 
M/R = 0.44, approaching the Buchdahl limit M/R -> 4/9. 
The behavior of B and Cl is qualitative very different from 
the previous two cases. However, small-scale structure still 
appears at late times. 





FIG. 11: Snapshots of B (lower panel) and fi (upper panel) 
along the zu/R = 0.5 line for a relativistic star with M/R — 
0.44. Analysis with FFT of the time series reveals that the 
initial data consists of mainly three modes with frequencies 
fi = 0.0356, f 2 = 0.0622 and / 3 = 0.0979. Curves are plotted 
at time intervals At — 5 from t = to t = 140 (about 5 
oscillation periods of the fundamental mode) . 



B^ is increased by a factor of l/a for a given initial angu- 
lar velocity profile fi(0; zu, z) according to Eq. (150). Fol- 
lowing this analysis, we expect that for M/R = 0.3, the 
maximum value of B on the w = 0.5 line should be about 
1.9 times larger than the Newtonian value, and the oscil- 
lation period to be longer by about a factor of 2.5. Our 
numerical result gives a factor of 2 for the amplitude and 
2.5 for the oscillation period on the cylinder w/R = 0.5. 
This indicates that the background spacctimc metric for 
this star changes slowly along the poloidal field lines and 
so the behavior of magnetic braking roughly agrees with 
the above simple analysis. 

To explore strong gravity regime, we consider an ex- 
treme relativistic star with M/R = 0.44, approaching the 
Buchdahl limit M/R -> 4/9. Recall that at the Buchdahl 
limit, the pressure at the center of the star (r = 0) be- 
comes infinite, which causes the lapse function a vanish 
at the center. We therefore expect to see a very differ- 
ent evolution pattern of B and fi. Figure 10 shows the 
snapshots of B and fi along the 9 = 45° line. We see 
that although the pattern of the curves is very differ- 
ent from the previous two cases, small-scale structures 
still develop at late times. The amplitude of the toroidal 
field is much larger and the oscillation timescale is much 
longer at small r than that at large r. It follows from 
Eqs. (130) and (131) that we can define the effective lo- 
cal Alfven speed as 



v e s(r) = 




(186) 



1 - 




FIG. 12: Same as Fig. 6 but for a relativistic star with M/R = 
0.44. 



Hence the Alfven speed is very small at the center when 
the Buchdahl limit is approached, which explains why 
the oscillation period is long close to the center. 

Figure 11 shows the snapshots of B and fi on the 
zu/R = 0.5 plane. For this extreme relativistic star, the 
solution of the MHD equations with the initial angular 
velocity distribution (172) is no longer close to the fun- 
damental mode. We perform a fast Fourier transform 
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(FFT) on the time series B(t;0.5,u) for three different 
values of u. We find that each of the three sets of time 
series consist mainly of three frequencies /i = 0.0356, 
f 2 = 0.0622 and f 3 = 0.0979. We also see that the 
peak of the toroidal field shifts with increasing time to a 
lower value of z. This can be explained by the following 
physical argument. The initial angular velocity can be 
regarded as a superposition of waves of B^ traveling to 
z = —zr = — V-R 2 — ra? 2 and z = zr. The waves are 
reflected at the two points and travel backward. Since 
the value of r is smaller at smaller z, it follows from 
Eq. (186) that the local Alfven speed decreases as z de- 
creases. However, the equatorial symmetry condition on 
O forces the toroidal field to vanish at the equator. As 
a result, the waves pile up near the equator, causing the 
peak to shift to a smaller value of z. 

From Eq. (180) we found O mea n = 0.64. Even though 
the ft in Fig. 11 consists of several modes, it still oscillates 
about this mean angular velocity. As in the previous 
cases, the rotational kinetic energy associated with the 
differential rotation transfers back and fro to the energy 
associated with the toroidal magnetic field (see Fig. 12). 



2. 1 = 2 initial field (Cases Ila-IIc) 

For the I = 2 field, we numerically integrate Eqs. (164) 
and (165) in the (u,v) coordinate system introduced in 
Section VB. The advantage of using these rather com- 
plicated coordinates is twofold. First, we only need to 
integrate a set of decoupled 1+1 evolution equations. 
Second, more importantly, small-scale structures will de- 
velop at late times in a direction across the poloidal 
field lines. Using other coordinate systems (e.g. stan- 
dard Cartesian, spherical or cylindrical coordinates) will 
involve taking derivatives across the poloidal field lines, 
which will be numerically inaccurate when the finite-sized 
grid can no longer resolve the growth of small-scale struc- 
tures. We will demonstrate this numerical inaccuracy in 
Section VI B. To estimate the numerical error of our 
code, we monitor the conserved integrals t2(v) and jziv) 
defined in Eqs. (168) and (171) as well as second order 
convergence. We find that the fractional variation of e 2 
and j2 is less than 0.3% for all of our runs with our canon- 
ical grid allocation. 

In the literature, the rotation law 



= A 2 (n c - O) 



(187) 



is often used to model differentially rotating relativistic 
stars (see, e.g. Ref. [62]). Here A is a constant. This law 
corresponds to, using the metric (89), the relation 



1 



1 + r - 



a 2 6 



(188) 



2 A 2 



where A = A/R and we have ignored the frame dragging 
term for simplicity. In the Newtonian limit (a — ► 1), Q p 



reduces to the so-called "j-constant" law [63]. Unfortu- 
nately, this rotation law does not satisfy the constraint 
B^(0)djfl = at r = R. To "fix" this, we choose a 
modified rotation profile [71] 

O(0; r, 9) = Cl p (f, 6)f(r) + 6 (f, 9)[1 - /(f)] , (189) 

where fi^ is a function that satisfies the constraint at r — 
R and / is a function which is close to 1 and drops rapidly 
to at f = 1. We also require /'(l) = 0. Specifically, we 
choose 

/(') = ^#7^^ d90) 



M°) - /f(1) + |/f(1) 



Mr) = 



1 + exp[(f - b)/e] 



(191) 



where b and e are constant parameters which we set to 
be b = 0.8 and 1/e = 25. We set 



n b (r,0) = (fl p (l,O))f(2-r) 



(192) 



where 



(1,9)) = — d4> d9sm9(l p (l, 9) (193) 
4?r J J 



(Cl } 



= § - tanh 1 



(194) 



We set the parameter A = 1 in all of our numerical cal- 
culations. The value of 0(0; r, 9) chosen in this way is 
close to O p (r, 9) except near the surface of the star, where 
the modification makes it satisfy the boundary condition 
B^(0)djQ = at r = R. 

Figure 13 shows the evolution of the toroidal field and 
angular velocity along the line w/R = 0.5 for a New- 
tonian star (M/R <C 1). Magnetic braking causes the 
toroidal field and angular velocity to oscillate. We see 
that small-scale structures develop at late times, which 
is expected since the direction does not line up with 
any poloidal field line, while coherent oscillations are re- 
stricted to field lines. 

Figure 15 shows the snapshots of B and O along the 
line v = 0.1 (sec Fig. 14) for the same star. Small-scale 
structures are not seen along this line since we are looking 
in the direction of a poloidal field line. We perform FFT 
on the time series of three chosen points on the line. We 
find that the FFT spectra consist of the same frequencies 
for the three points. The lowest four cigenfrequencies 
are determined to be 0.35, 0.65, 0.95 and 1.25 in our 
nondimcnsional units. 

Similar to the I = 1 cases, we can define, from 
Eq. (171), a mean angular velocity at each v by 



O mC an(£) = h{v)/I 2 {v) , 



where 



h{v) = / du 



sin 



o a 2 A(3cos 2 6> + 1) 



(195) 
(196) 
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FIG. 13: Snapshots of B (lower curves) and ti. (upper curve) 
along the w/R = 0.5 line for a Newtonian star (M/R <C 1). 
The initial magnetic field is the I = 2 field given by Eq. (117). 

The surface of that star is located at z = z R = R? - vj 2 . FIG - 14: The v = 

Since the direction does not line up with the field lines, the 
build-up of small-scale structure is observed, as expected. 



For this Newtonian star, we found fi mcan (0.1) = 0.70. We 
see from Fig. 15 that Cl oscillates about this mean value. 
Figure 16 shows the evolution of the reduced energies 
at v = 0.1. As expected, we see that both the reduced 
rotational kinetic energy e™* and magnetic energy e™ ag 
oscillate, while keeping the total energy i 2 conserved. 

Figure 17 shows the snapshots of B and ft along the 
zu/R — 0.5 line for a relativistic star with M/R = 0.3. 
The graphs look qualitatively the same as the Newto- 
nian situation. Similar to the I = 1 case, it is easy to 
see from Eqs. (164) and (165) that if the relativistic fac- 
tors a and A do not change significantly inside the star, 
the main relativistic effects are to cause time to dilate 
according to Eq. (185) and to increase the amplitude of 
the toroidal field by a factor of 1/a for a given set of 
initial data. This explains the similarity of the features 
in Figs. 13 and 17, and that small-scale structures ap- 
pear later for the relativistic star. Figure 18 shows the 
snapshots along the v = 0.1 line, and Fig. 19 shows the 
evolution of the reduced energies. The angular veloc- 
ity function Ct(t;u, 0.1) oscillates about the mean value 
^mean(O.l) = 0.45. The FFT spectrum reveals the lowest 
four eigenfrcquencies to be 0.14, 0.26, 0.38 and 0.62. Note 
that these frequencies are 40% of those of the Newtonian 
stars, which agrees with the time dilation effect predicted 
by Eq. (185). However, the ratio of the amplitudes of the 
toroidal field is much larger than the predicted value a. 
This is because the two stars do not have the same initial 
angular velocity profile [see Eqs. (188) and (189)]. 

Finally, we explore effect of strong gravity by consider- 
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FIG. 15: Snapshots of B (lower panel) and ft (upper panel) 
along the v = 0.1 line for a Newtonian star (M/R <C 1). The 
FFT spectrum indicates that the first four eigenfrequencies 
to be approximately 0.35, 0.65, 0.95 and 1.25. Curves are 
plotted at time intervals At = 0.5 from t = to t = 15 
(about 5 oscillation periods of the fundamental mode). 



ing a relativistic star with M/R = 0.44, approaching the 
Buchdahl limit. Figures 20 (21) shows the snapshots of 
B and (l along the zu/R = 0.5 (v = 0.1) line. Figure 22 
shows the evolution of the reduced energies at v = 0.1. 
We see that the behavior is qualitatively very different 
from the previous two stars. Small-scale structures still 
develop, but it happens at much later time. The oscil- 
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FIG. 16: Evolution of the reduced energies for a Newtonian 
star with I = 2 initial magnetic field. The solid line is e!j ot (i) — 
£2 Ct (0), the dotted line is e™ ag (i), and the dashed line is £2 — 
e 2 ot (0). All the reduced energies are evaluated at v = 0.1. 




U 



FIG. 18: Snapshots of B (lower panel) and Q, (upper panel) 
along the v — 0.1 line for a relativistic star with M/R — 0.3. 
The FFT spectrum reveals that the first four eigenfrequencies 
to be 0.14, 0.26, 0.38 and 0.62. Curves are plotted at time 
intervals At = 1 from t = to t = 35 (about 5 oscillation 
periods of the fundamental mode). 
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FIG. 17: Same as Fig. 13 but with M/R = 0.3. 



lation amplitudes of B and Cl are much larger along the 
{inconstant lines. 

The growth of the toroidal field twists the poloidal field 
lines in the azimuthal direction. The magnetic field lines 
are therefore still confined to the surface spanned by the 
poloidal field lines. However, since we consider the sit- 
uation in which the rotational timescale is much longer 
than the Alfven timescale, the induced toroidal field is 
very small compared to the poloidal field. In order to see 
the small twist of the field lines, we multiply the toroidal 
field by an arbitrary factor of 10va/Q c R in Fig. 23, 
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FIG. 19: Same as Fig. 16 but for a relativistic star with 
M/R = 0.3. 



where the field lines on the v = 0.1 surface is plotted 
for a relativistic star with M/R — 0.3. As expected, the 
oscillations of the toroidal field lines cause the total field 
lines to twist back and forth in the azimuthal direction. 
The field lines of relativistic stars with other compaction 
factors M/R behave in a similar way. 
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FIG. 20: Same as Fig. 13 but with M/R = 0.44 (The upper 
curve is fi + 2). 
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FIG. 22: Same as Fig. 16 but for a relativistic star with 
M/R = 0.44. 
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FIG. 21: Snapshots of B (lower panel) and Cl (upper panel) 
along the v — 0.1 line for a relativistic star with M/R = 0.44. 
The FFT spectrum reveals that the first four eigenfrequencies 
to be 0.016, 0.031, 0.045 and 0.060. Curves are plotted at time 
intervals At = 15 from t = to t = 300 (about 5 oscillation 
periods of the fundamental mode). 
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B. Role of Viscosity (Case III) 

Our previous results show that the phase mixing aris- 
ing from magnetic braking is likely to induce irregular 
angular velocity flow. A significant meridional current 
will build up eventually, and thus may induce other mag- 
netic instabilities and cause turbulence. Magnetic diffu- 
sion and viscosity will play a significant role on the subse- 



FIG. 23: The magnetic field lines as functions of time along 
the axisymmetric surface v = 0.1 for a relativistic star with 
M/R = 0.3. The poloidal field lines are the same as the field 
lines at t = 0. The toroidal field has been multiplied by an 
arbitrary factor of IOva/^cR so that they are more visible. 
The circles on the upper and lower ends are the boundary of 
the star, i.e., the lines of revolution generated by the points 
ui and U2 in Fig. 14. 
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FIG. 24: Energy and angular momentum versus time by 
evolving Eqs. (146) and (147) for the I = 2 poloidal field 
using the cylindrical-like coordinates (uj, u) in the absence of 
viscosity. The compaction M/R is 0.3. The results for two 
resolutions (32 x 32 and 64 x 64) are plotted. The energy and 
angular momentum derivate from their initial values at late 
times because the finite-sized grid fails to resolve the small- 
scale structures built up by magnetic braking. Increasing the 
resolution postpones, but does not eliminate, the breakdown 
of conservation. 

quent evolution. It has been shown in Newtonian MHD 
that magnetic diffusion will damp out the angular veloc- 
ity oscillations along a poloidal line. The final state of 
the star will have a constant angular velocity profile on 
each magnetic surface (see [67] and references therein). 
Here we wish to discuss the role of viscosity in damp- 
ing differential motion. Microscopic viscosity in a typical 
neutron star is very small, and its timescale is very long 
compared to the rotational period [64] . Here we may re- 
gard viscosity as a "turbulent viscosity," which models 
the effects of turbulence via an effective shear viscosity. 
The typical turbulent viscosity is v ~ lAv (see Ref. [55], 
§ 33), where I is the length scale of the largest turbu- 
lent eddies and Av is the velocity fluctuation over the 
distance I. This turbulent viscosity is much larger than 
the microscopic viscosity, and can affect the flow on a 
dynamical timescale. 

We solve the MHD equations (164) and (165) with vis- 
cosity v ^ 0. We only present one case here, as the result 
is qualitatively similar for all the other cases. The initial 
poloidal field is given by the I = 2 field and the M/R 
of the star is 0.3. The viscosity coefficient v is set to be 
0.002. This value is chosen so that the toroidal field and 
angular velocity can go through several Alfven oscilla- 
tions before the viscosity damps the oscillations. 

As mentioned in Section VB, the (u,v) coordinates 
are not convenient when viscosity is present. We use 
the cylindrical- like coordinates (w,u) introduced in Sec- 



FIG. 25: Same as Fig. 24 but with viscosity v = 0.002. Viscos- 
ity suppresses the formation of small-scale structures. Hence 
the integrals can be conserved accurately even in cylindrical 
coordinates. 

tion VA in this case. In the absence of viscosity, we are 
not able to integrate the MHD equations accurately in 
these coordinates very long because our finite-sized grid 
can no longer resolve the small-scale structures that de- 
velop at late times. Figure 24 illustrates the resulting 
numerical inaccuracy by plotting the energy E and an- 
gular momentum J computed from the numerical data as 
a function of time. We evolve the MHD equations with 
resolutions 32 x 32 and 64 x 64. We see that the two 
conserved integrals derivate from their initial values at 
late times. Doubling the resolution does not eliminate, 
but only postpones, the breakdown of E and J conserva- 
tion. Hence in the absence of viscosity, we must use the 
(u, v) coordinates to ensure a stable evolution. However, 
when we include a small viscosity v = 0.002, the situ- 
ation changes drastically. The viscosity suppresses the 
build up of small-scale structure. As a result, the numer- 
ical inaccuracy disappears (see Fig. 25). The fluctuation 
of energy E is 3% (0.4%) in 32 x 32 (64 x 64) resolution. 
The fluctuation of angular momentum J is 1% (0.4%) in 
32 x 32 (64 x 64) resolution. 

Figure 26 shows the snapshots of B and ft at the equa- 
tor. We see that the star is driven to uniform rotation 
with no toroidal field at late times, as expected. Since 
angular momentum is conserved, we can calculate the fi- 
nal angular velocity from Eq. (153), which is also valid 
for the I — 2 field. We obtain 

Ofi na] = J/J, (197) 

f f du 

I = 3a R / dww 3 y/l - w 2 \ . (198) 

Jo Jo a * 

Not surprisingly, our numerical value of rifi na i shown in 
Fig. 26 agrees with the above formula, f2fi na i = 0.425. 
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FIG. 26: Snapshots of B (solid lines) and Q (dashed lines) 
at the equator in the presence of viscosity. The coefficient of 
viscosity is set to be v — 0.002. The toroidal field vanishes and 
the star becomes rigidly rotating at late times, as expected. 
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FIG. 27: Evolution of energies for a relativistic star with 
M/R = 0.3 in the presence of viscosity. The solid line is 



E Iot (t) 



E ro t(0), the dotted line is .E mag (t), the dot-dashed 



line is E vla (t), and the dashed line is the sum of the three 
lines, i.e. E rot (i) + £ mag (t) + E vis (i) - E rot (0), which should 
remain constant and equal to 0. The deviation from zero of 
the dashed curve is due to truncation error of our code, which 
decreases with increasing resolution. 



Finally, Fig. 27 shows the evolution of E rot (i), E mag (i), 
and E v i s (t). We see that a portion of rotational energy 
transfers to magnetic energy of the toroidal field. Both 
energies are eventually turned into heat. The total energy 
is conserved up to the truncation error, which decreases 
with increasing resolution. 



VII. SUMMARY AND DISCUSSION 

We have studied the magnetic braking in differen- 
tially rotating, relativistic stars. The results we found 
are very different from the models studied in Papers I 
and II. The models in Paper I are incompressible and one- 
dimensional (infinite cylinders), the models in Paper II 
are compressible and one-dimensional (infinite cylinders), 
while the models in this paper are incompressible and 
two-dimensional (axisymmetric, finite stars). In Papers I 
and II, the Newtonian MHD equations were evolved while 
here we evolved the relativistic MHD equations with the 
assumption that the initial poloidal fields is large com- 
pared to rotational energy, but weak compared to the 
gravitational binding energy. In Papers I and II, it was 
found that magnetic braking generates coherent toroidal 
magnetic fields and laminar velocity flows. In the mod- 
els considered here, we found that the field breaks down 
because of the phase mixing both in Newtonian and rel- 
ativistic regimes. Our analysis indicates that the phase 
mixing is a result of the geometry of the stars and the 
poloidal field lines: phase mixing is absent in the infinite 
cylinders of Papers I and II simply because the systems 



are highly symmetric. However, even in finite stars, if 
viscosity is sufficiently strong, the toroidal fields may be 
damped before phase mixing becomes appreciable. 

Differential rotation is damped by the presence of 
poloidal magnetic fields and viscosity. Magnetic brak- 
ing causes the toroidal fields and angular velocities to 
undergo oscillations along each poloidal field line. In the 
regime explored in this paper, the oscillations along a 
given field line are independent of those along other lines. 
This phase mixing effect is likely to stir up turbulent-like 
motion. When viscosity is present, the star ultimately 
is driven toward uniform rotation. When the compact- 
ness of the star, M/R, is not very large, the effect of 
magnetic braking in a relativistic star is similar to the 
Newtonian case, but the timescale of the braking process 
is increased roughly by a factor (1 — 2M/R) 2 . When M/R 
approaches the Buchdahl limit 4/9, the braking process 
is strongly affected by the spacetime curvature. Strong 
toroidal fields pile up at small r along each field line, 
since the proper time elapses more slowly than at larger 
r for the same amount of coordinate time. The whole 
process of magnetic braking can take much longer in this 
case than in the Newtonian situation. 

The assumptions necessary to bring the 2+1 MHD 
equations into a set of 1+1 equations are that the sys- 
tem is axisymmetric; viscosity, magnetic diffusion and 
meridional currents can be neglected; and both magnetic 
and rotational energies are much smaller than the grav- 
itational binding energy. The condition E mag » E rot 
is a sufficient but may not be a necessary condition for 
the meridional currents to remain small in tens of Alfven 
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times. It is possible that the meridional currents can still 
be small compared to the rotational velocity even if this 
condition is not satisfied. If this is the case, our analysis 
here will still be valid. It is possible that most young 
neutron stars or merged binary neutron star remnants 
do not have strong poloidal fields and the opposite limit 
E Iot S> E mag may be more relevant. More detailed MHD 
calculations are required to study the general case. 

The microscopic viscous timescale is of order t v = 
R 2 /Sv w 10 9 s for a typical neutron star, where v = 
cm 2 s 1 and T is temperature [64]. How- 
ever, turbulent viscosity can be much bigger and might 
drive the star toward uniform rotation on much shorter 
timescale. The Alfven timescale in a neutron star with 
magnetic field B = 10 12 G is tens of seconds (Paper I). 
Hence, the effect of magnetic field is much more impor- 
tant than microscopic viscosity in neutron stars. Mag- 
netic braking may have important consequences for gravi- 
tational wave signals and gamma-ray bursts, as discussed 
in the Introduction. 

A significant amount of meridional current may also 
be generated as a result of the phase mixing. Our cal- 
culations do not take into account the later build-up of 
meridional currents, which may induce convective insta- 
bilities [65], driving the seed magnetic field to high values 
greatly exceeding 10 12 G [66]. The meridional currents 
may also induce other possible MHD instabilities, which 
may contribute to the redistribution of angular momen- 
tum [53, 67]. Our goal here is to show that even in the 
simplest case, magnetic braking can also induce irregular, 
turbulent-like behavior. 

More realistic evolutionary calculations of magnetic 
braking in neutron stars should clarify some of the above 
issues. One computational subtlety is that the Alfven 
timescale is usually much longer than the dynamical 
timescale of the star. In this regime it may prove too 
taxing to a relativistic MHD code to evolve a differen- 
tially rotating star for the required length of time for 
magnetic braking to take effect. One possibility is to 
treat part of the evolution in the quasistatic approxi- 
mation, as in a typical stellar evolution code, up to the 
moment that stable equilibrium can no longer be sus- 
tained. One other possibility is to artificially amplify the 
magnetic field so that the effect of magnetic braking will 
show up in a computationally managable timescale, and 
then scale the results for smaller ratios. Still another ap- 
proach is to use implicit differencing scheme to avoid the 
Courant constraint on the evolution timestep. However, 
our calculations suggest that small-scale irregular angular 
velocity flows and meridional currents are likely to grow 
during the process of magnetic braking. The ability of 
a numerical code (finite difference or spectral) to resolve 
this behavior on ever decreasing scales is an important 
challenge that will require further analysis. We hope to 
tackle the more general problems in the near future by 
means of numerical simulations in full general relativity. 
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APPENDIX A: VALIDITY OF IGNORING v r 

AND v e 



The main effect of magnetic braking is to wind up 
the frozcn-in magnetic field and generate toroidal fields, 
which back- react on the angular velocity, = O. The 
change in angular velocity ultimately drives a meridional 
current. In this appendix, we estimate the timescale for 
v r and v e to become comparable to ?A Here we only 
study the simplest case: a differentially rotating, incom- 
pressible Newtonian star, with an initial I — 1 poloidal 
field .6(0) = B n e Z} and zero viscosity. 

We start from the Newtonian MHD equations (1) 
(4). We neglect viscosity in this appendix. We find it 
convenient to work in cylindrical coordinates (m,$,z). 
The meridional components of the velocity are w ro = 
sin0v f + cos6v®, and v z = cos9v f - sin0v § . The MHD 
equations become 



d t B* 

m 

d t B* 



B d z n + {(b z - B )d z n + B^d^n 

-v>djB*] , (Al) 



4irp 



— ( B^d^B^ + —B^B^ 
Airp \ w 



-v 3 dM v^VL 

w 

= I)' (),!■' - r* {/J!' , 
= -v 3 d.v™ + vjVl 2 



(A2) 
(A3) 



-d m p-d m $ 



p 



-J-a^ + ^d^-^r] (A4) 

d t v z = -v j div z - -d z P - d z <P - —d z B 2 
p 8irp 

+ ^ 9j B z , (AS) 



d t B™ 
d t B z 



v J djB™ - B 3 djV u 
v j djB z - B 3 djv z 



(A6) 
(A7) 



where i and j denotes zu and z, and the summation con- 
vention is assumed. At t = 0, we assume w ro (0;ra7, z) = 
v z (0;w,z) = 0, ^(0;n7,z) = Q(0;w,z), B(0;w,z) = 
Boe z , and the star is in hydrostatic equilibrium: 

mn 2 (0;m, z) - i<9 ro P(0;n7, z) - 9 ro $(0;n7, z) = 0, (A8) 

--d z P(0;w,z)-d z $(0;zn,z) = 0, (A9) 
P 

where we have used Eqs. (A4) and (A5). After grouping 
these equilibrium terms, we see that the only remaining 
nonvanishing term at t — in Eqs. (A1)-(A7) is B d z tt, 
the first term on the right side of Eq. (Al). This simply 
tells us the obvious fact that when there is differential 
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rotation along the initial field lines (which is in the z- 
direction for this initial field), a toroidal field will be 
generated by a term linear in 0. It follows from Eq. (A2) 
that the toroidal field will then modify O. Given our as- 
sumptions of slow rotation and weak magnetic field, all 
the other terms are of higher order, at least at small t. 
This is the reason we neglect them in this paper. How- 
ever, it is possible that u ro and v z will become important 
at late times. 

Here we crudely estimate the timescale in which w ro 
and v z are driven to a magnitude comparable to by 
the change of induced by magnetic braking. We fo- 
cus on Eq. (A4). We assume that the pressure P and 
gravitational potential <I> in our nearly spherical configu- 
ration do not change and set w ro = v z = 0, P ro = 0, and 
B z = B on the right side of Eq. (A4). Using Eq. (A8), 
Eq. (A4) becomes 



time. This crude calculation therefore does not take into 
account the possibility of the growth of other instabilities, 
like the magnetorotational instability (MRI). However, 
the growth time of this instability is also of order P rot (sec 
e.g., Ref. [53], Section IV). Hence we conclude that our 
numerical results are valid as long as the Alfven timescale 
tA < -Prot, or equivalently P rot < Pmag- 

Although our analysis in this paper is restricted to the 
early phases of the evolution P rot 3> t ^ tA in the slow 
rotation, weak magnetic field limit, we are able to deter- 
mine the nonlinear evolution of the angular velocity and 
toroidal magnetic field profiles during this period for a 
fully relativistic, differentially rotating star. 



APPENDIX B: TRANSFORMATION BETWEEN 
u, v, r AND FOR I = 2 FIELD 



d t v™ = vj[n 2 (t) - o 2 (o)] - . (aio) 

We then substitute for O and B^ the analytic solutions 
given by Eqs. (55) and (56). For simplicity, we assume 
that only one mode is present, and write 



n(t) = O(0) + e ^1 - J cos k n z (1 - cosLO n t) , (All) 

Bf(t) = B (jj^j (l - ^ sink n z sinw„i , (A12) 

where e is a constant parameter of dimension 1/time, 
measuring the degree of differential rotation along the 
initial magnetic field lines. We assume that the star may 
have a high degree of differential rotation so that e can 
be of the same order as 17(0). The factor (1 — w 2 /R 2 ) is 
inserted to make both O and B^ regular at w = R. Equa- 
tion (A10) can now be solved analytically. The growth 
of w ro is dominated by the terms that go linearly in time: 



tzue 
+2 [ I 



R 2 
R? 



2O(0) 



cos 2 k„ z 



cos k n z 
1 

" 2 



P 2 " 



. (A13) 



Hence w ro will be comparable to ~ fl(0)w when 
t > t c ~ O(0)/e 2 > 1/O(0). Hence t c > P Iot , where 
P ro t is a rotation period. It follows from the incompress- 
ibility condition V • v = that v z is in general of the 
same order as . Therefore, v z will also develop on 
the timescale t c . On the other hand, the timescale for 
magnetic braking to affect and generate B^ is tA, the 
Alfven timescale. Although our calculation is done on a 
Newtonian model with an I = 1 initial field, we expect 
that the result is similar in general relativity and with 
other initial magnetic field configurations. 

Our estimate is based on dropping all the terms on the 
right side of Eq. (A4) which are of higher order at early 



To solve Eqs. (164) and (165), we need to know the 
transformation between r, 9 and u, v. In this appendix, 
we show the required transformation. 

It is convenient to introduce the cylindrical variables 

w = fsin9 , z — f cos 8, (Bl) 
where f — r/R. We also define an auxiliary variable 

= u/R 2 = -f 2 (3 cos 2 6>- 1) . (B2) 
The relations between w, z, and u*, v are 



137 — ZZ 



3V3 



w 2 z . 



(B3) 
(B4) 



For a given value of v, we denote by w\ and ra7 2 (with 
V02 > mi) the two points at which the v =constant line 
intercepts the surface of the star f = 1. Substituting 
z = VT — vo 2 in Eq. (B4), we obtain the equation for wi 
and W2- 



v = 



3\/3 



-w 



(B5) 



The equation can be turned into a cubic equation and 
the solution is given by 



VJi 




'-(1 + 2 008^) 



-sin v . 
3 



(B6) 

(B7) 
(B8) 



It follows from Eqs. (B3) and (B4) that the corresponding 
Mi* and «2* are given by 



8v 2 
27^ 



, (.3 = 1, 2) . 



(B9) 



25 



The transformation from (r, 9) to (u, v) can be per- 
formed by the following procedure. We first compute 
and v from Eqs. (B1)-(B4). We then compute uu and 
u 2 * from Eqs. (B6)-(B9). Finally, we compute u by the 
formula [see Eq. (162)] 



u = 



11* Wl* 
«2* — Uu 



(BIO) 



To transform from (u, v) to (r, 9), we first compute Ui* 
and w 2 * from Eqs. (B6)-(B9). Next we calculate from 
Eq. (BIO). We then solve Eqs. (B3) and (B4) to obtain w 
and z. The two equations can be combined into a cubic 
equation, which is solved by the standard cubic equation 
formula. Finally, the values of r and 9 is obtained from 
Eq. (Bl). 



APPENDIX C: FINITE DIFFERENCING 
SCHEME 

In this appendix, we describe in details our finite dif- 
ferencing scheme for the I = 1 MHD equations (151) 
and (152) in the presence of viscosity. The other cases 
are treated in a similar fashion. 

We define 

X™. = X(t n ;u i ,w j ) , (CI) 

where X is any variables and 

t n = nAi , u t = iAu , wj = jAw . (C2) 

We use a uniform spatial grid (Aw = constant, Am — 
constant) and a uniform time step (At =constant) in our 
calculations. We store B at grid points (i,j) and Cl at (i+ 
1/2, j). We compute the right side of Eqs. (151) and (152) 
by the following finite difference approximation: 

d i B(f l ;u l ,m J ) « f B ;ij{& n ) , (C3) 
dp,(t n \u l+l/ 2,i5jj) w /n ; i+i/2,j(-B") 

+U.,+i/2.jm > (C4) 



where 

fn-i+i/2,j(B n ) = 



\ / O™ — o™ 

a i,j A i,3 / "*+l/2,j "t-l/2,j 



^7 



Aw 



Aw 



(C5) 



|(C6) 



We use the expression in Eq. (149) to com- 
pute the finite difference term, f V ;i+i/2,j(& n ), f° r 
Ui+1/21 ^')]vis- Specifically, we first calculate 
d„f2, c^fl, c?^Sl and d^d-ati by standard central 

differencing, and then compute dffl, dffl, dgCl and d 2 Cl 
by the transformation formulae 



r r(l — w 1 ) 



(C7) 



C^f) = M\/ 1 — tI7 2 9ro^ — 



\/l — 1Z7 ; 



:(l-{t 2 )^^ (C8) 



d 2 n 



3zu 2 u 
+ f 2 (l-w 2 ) 2 

-.2/1 ~2\2a2 



2a7U 



f 2 (l-II7 2 ) 2 " f 2 (l-tU 2 ) 



-d A d f ,fl 



(C9) 



9 2 fi = ^(l-^)^0 + T ^(l-^)^0 
—2mu(l — u 2 )dtidu£l — wd^Cl 

..112m 2 , a 
-w(l - u )- r^-«fi^ . 



(CIO) 



Finally, we compute f V ;i+i/2,j(& n ) using Eq.. (149). 

We use an iterated Crank-Nicholson [57] scheme to 
evolve Eqs. (C3) and (C4). To be specific, we first predict 
B and f2 at the next time step by 

^B^ 1 = B^ + Atf B .ij(h n ) , (Cll) 

(1) ^ + + i/ 2 , = ^ + i/2j+At[/ Bii+1/2J (B n ) 

+/ Wi i+i/2j(A n )] . (C12) 

Next we estimate B and f2 at the next half time step by 



(l)An+l/2 1 

' ~" 2 

(l)6«+l/2 _ 1 

"i+l/2,j - 2 



R™ -4- W R™ +1 



"i+l/2,j T t+l/2,j 



(C13) 
(C14) 



We then use these variables to obtain a more accurate 
estimate of B and fl at the next time step: 

(2)^+! = B? d + Atf B ;ij( {1) & n+1/2 ) , (C15) 

^ +1/2j +At[/B; J+ l/ 2j ( (1) S" +1 / 2 ) + 

U, + i/2,( {1) n n+1/2 )} ■ (ci6) 



(2)^n 



+ 1 

+ 1/2J 



To ensure the stability of this finite differencing scheme, 
we have to do this "corrector" step one more time [57]. 
Hence we compute 



(2)An+l/2 _ 1 

(2)A«+l/2 1 
"i+l/2,j - 2 



R» i ( 2 ) R™+! 



6™ , +( 2 )6" +1 

il i+l/2,j + S S+l/2j 



(C17) 
, (C18) 



and use them to obtain the final values of B and f2 at 
the next time step: 

Spt 1 = +Ai/ B;JJ (( 2 )r! n + 1 / 2 ) , (C19) 

A ?+1/2J = ^ + l/2,,+At[/ B;i+1/2j (( 2 ^" +1 / 2 ) 

+ /,;,+l/2 J ( (2) ^ n+1/2 )] • (C20) 

The iterated Crank-Nicholson scheme is second order ac- 
curate in space and time. 
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